home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_K1.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.8 KB  |  221 lines
  1. /* specfunc/bessel_K1.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_bessel.h>
  27.  
  28. #include "error.h"
  29.  
  30. #include "chebyshev.h"
  31. #include "cheb_eval.c"
  32.  
  33. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  34.  
  35. /* based on SLATEC besk1(), besk1e() */
  36.  
  37. /* chebyshev expansions 
  38.  
  39.  series for bk1        on the interval  0.        to  4.00000d+00
  40.                     with weighted error   7.02e-18
  41.                      log weighted error  17.15
  42.                    significant figures required  16.73
  43.                     decimal places required  17.67
  44.  
  45.  series for ak1        on the interval  1.25000d-01 to  5.00000d-01
  46.                     with weighted error   6.06e-17
  47.                      log weighted error  16.22
  48.                    significant figures required  15.41
  49.                     decimal places required  16.83
  50.  
  51.  series for ak12       on the interval  0.        to  1.25000d-01
  52.                     with weighted error   2.58e-17
  53.                      log weighted error  16.59
  54.                    significant figures required  15.22
  55.                     decimal places required  17.16
  56. */
  57.  
  58. static double bk1_data[11] = {
  59.    0.0253002273389477705,
  60.   -0.3531559607765448760, 
  61.   -0.1226111808226571480, 
  62.   -0.0069757238596398643,
  63.   -0.0001730288957513052,
  64.   -0.0000024334061415659,
  65.   -0.0000000221338763073,
  66.   -0.0000000001411488392,
  67.   -0.0000000000006666901,
  68.   -0.0000000000000024274,
  69.   -0.0000000000000000070
  70. };
  71.  
  72. static cheb_series bk1_cs = {
  73.   bk1_data,
  74.   10,
  75.   -1, 1,
  76.   8
  77. };
  78.  
  79. static double ak1_data[17] = {
  80.    0.27443134069738830, 
  81.    0.07571989953199368,
  82.   -0.00144105155647540,
  83.    0.00006650116955125,
  84.   -0.00000436998470952,
  85.    0.00000035402774997,
  86.   -0.00000003311163779,
  87.    0.00000000344597758,
  88.   -0.00000000038989323,
  89.    0.00000000004720819,
  90.   -0.00000000000604783,
  91.    0.00000000000081284,
  92.   -0.00000000000011386,
  93.    0.00000000000001654,
  94.   -0.00000000000000248,
  95.    0.00000000000000038,
  96.   -0.00000000000000006
  97. };
  98. static cheb_series ak1_cs = {
  99.   ak1_data,
  100.   16,
  101.   -1, 1,
  102.   9
  103. };
  104.  
  105. static double ak12_data[14] = {
  106.    0.06379308343739001,
  107.    0.02832887813049721,
  108.   -0.00024753706739052,
  109.    0.00000577197245160,
  110.   -0.00000020689392195,
  111.    0.00000000973998344,
  112.   -0.00000000055853361,
  113.    0.00000000003732996,
  114.   -0.00000000000282505,
  115.    0.00000000000023720,
  116.   -0.00000000000002176,
  117.    0.00000000000000215,
  118.   -0.00000000000000022,
  119.    0.00000000000000002
  120. };
  121. static cheb_series ak12_cs = {
  122.   ak12_data,
  123.   13,
  124.   -1, 1,
  125.   7
  126. };
  127.  
  128.  
  129. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  130.  
  131. int gsl_sf_bessel_K1_scaled_e(const double x, gsl_sf_result * result)
  132. {
  133.   /* CHECK_POINTER(result) */
  134.  
  135.   if(x <= 0.0) {
  136.     DOMAIN_ERROR(result);
  137.   }
  138.   else if(x < 2.0*GSL_DBL_MIN) {
  139.     OVERFLOW_ERROR(result);
  140.   }
  141.   else if(x <= 2.0) {
  142.     const double lx = log(x);
  143.     const double ex = exp(x);
  144.     int stat_I1;
  145.     gsl_sf_result I1;
  146.     gsl_sf_result c;
  147.     cheb_eval_e(&bk1_cs, 0.5*x*x-1.0, &c);
  148.     stat_I1 = gsl_sf_bessel_I1_e(x, &I1);
  149.     result->val  = ex * ((lx-M_LN2)*I1.val + (0.75 + c.val)/x);
  150.     result->err  = ex * (c.err/x + fabs(lx)*I1.err);
  151.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  152.     return stat_I1;
  153.   }
  154.   else if(x <= 8.0) {
  155.     const double sx = sqrt(x);
  156.     gsl_sf_result c;
  157.     cheb_eval_e(&ak1_cs, (16.0/x-5.0)/3.0, &c);
  158.     result->val  = (1.25 + c.val) / sx;
  159.     result->err  = c.err / sx;
  160.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  161.     return GSL_SUCCESS;
  162.   }
  163.   else {
  164.     const double sx = sqrt(x);
  165.     gsl_sf_result c;
  166.     cheb_eval_e(&ak12_cs, 16.0/x-1.0, &c);
  167.     result->val  = (1.25 + c.val) / sx;
  168.     result->err  = c.err / sx;
  169.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  170.     return GSL_SUCCESS;
  171.   }
  172. }
  173.  
  174.  
  175. int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result)
  176. {
  177.   /* CHECK_POINTER(result) */
  178.  
  179.   if(x <= 0.0) {
  180.     DOMAIN_ERROR(result);
  181.   }
  182.   else if(x < 2.0*GSL_DBL_MIN) {
  183.     OVERFLOW_ERROR(result);
  184.   }
  185.   else if(x <= 2.0) {
  186.     const double lx = log(x);
  187.     int stat_I1;
  188.     gsl_sf_result I1;
  189.     gsl_sf_result c;
  190.     cheb_eval_e(&bk1_cs, 0.5*x*x-1.0, &c);
  191.     stat_I1 = gsl_sf_bessel_I1_e(x, &I1);
  192.     result->val  = (lx-M_LN2)*I1.val + (0.75 + c.val)/x;
  193.     result->err  = c.err/x + fabs(lx)*I1.err;
  194.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  195.     return stat_I1;
  196.   }
  197.   else {
  198.     gsl_sf_result K1_scaled;
  199.     int stat_K1 = gsl_sf_bessel_K1_scaled_e(x, &K1_scaled);
  200.     int stat_e  = gsl_sf_exp_mult_err_e(-x, 0.0,
  201.                                            K1_scaled.val, K1_scaled.err,
  202.                        result);
  203.     result->err = fabs(result->val) * (GSL_DBL_EPSILON*fabs(x) + K1_scaled.err/K1_scaled.val);
  204.     return GSL_ERROR_SELECT_2(stat_e, stat_K1);
  205.   }
  206. }
  207.  
  208. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  209.  
  210. #include "eval.h"
  211.  
  212. double gsl_sf_bessel_K1_scaled(const double x)
  213. {
  214.   EVAL_RESULT(gsl_sf_bessel_K1_scaled_e(x, &result));
  215. }
  216.  
  217. double gsl_sf_bessel_K1(const double x)
  218. {
  219.   EVAL_RESULT(gsl_sf_bessel_K1_e(x, &result));
  220. }
  221.